%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Monte-carlo wavefunction simulation for a two-level atom [J. Dalibard et al., PRL68,0580 (1992)].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MonteCarlo
Constant;
figure(1);						clf;

G =2*pi*6e6;
W =     5*G;
D =    .5*G;
dt=   .5e-9;
t =[0:400]*dt;

M =MonteCarlo_Bloch(G,W,D);
rho0=[1 0 0 0]';
for jj=1:length(t)
	rho=expm(M*t(jj))*rho0;
	rh1(jj)=real(rho(1));
end

H =MonteCarlo_Hamilton(G,W,D);
nn1=zeros(size(t));
for jj=1:200
	Y =[1 0]';
	for j=1:length(t)
		n1(j)=Y(1)'*Y(1);  

		Y=expm(-i*H*dt)*Y;		%Y=(eye(2)-i*H*dt)*Y;
		Z=rand;
		if Y'*Y>Z
			Y=Y/(Y'*Y);
		else
			Y=[1 0]';
		end
	end
	nn1=nn1+n1;
   
	subplot(211),plot(t*1e9,rh1,'k',t*1e9,n1,'g','LineWidth',1.4)
	set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'FontSize',16,'LineWidth',1.4)
	ylabel('\fontsize{16}{\it\rho}_{11}')
	title(['Run = ',num2str(jj)]);
	subplot(212),plot(t*1e9,rh1,'k',t*1e9,nn1/jj,'b','LineWidth',1.4)
	set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'FontSize',16,'LineWidth',1.4)
	xlabel('\fontsize{16}{\itt}   (ns)');
	ylabel('\fontsize{16}{\it\rho}_{11}')
	pause(.01);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bloch equation matrix.

function M = MonteCarlo_Bloch(G,W,D)

M=[	0		G		i/2*W		-i/2*W
	0		-G		-i/2*W		i/2*W
	i/2*W	-i/2*W	-i*D-G/2	0
	-i/2*W	i/2*W	0			i*D-G/2	];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hamilton matrix.

function H = MonteCarlo_Hamilton(G,W,D)

H=[	0 	W/2
	W/2	-D-i*G/1	];
